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Abstract One of the most salient spatio-temporal patterns in population ecology is the synchroniza- 
pi I tion of fluctuating local populations across vast spatial extent. Synchronization of abundance has been 

D , widely observed across a range of spatial scales in relation to rate of dispersal among discrete populations. 

■ However, the dependence of synchrony on patterns of among-patch movement across heterogeneous land- 

scapes has been largely ignored. Here we consider the duration of movement between two predator-prey 
communities connected by weak dispersal, and its effect on population synchrony. More specifically, we 
introduce time delayed dispersal to incorporate the finite transmission time between discrete populations 
across a continuous landscape. Reducing the system to a phase model using weakly connected network 
theory, it is found that the time delay is an important factor determining the nature and stability of 
phase-locked states. Our analysis predicts enhanced convergence to stable synchronous fluctuations in 
general, and a decreased ability of systems to produce in-phase synchronization dynamics in the presence 
of delayed dispersal. These results introduce delayed dispersal as a tool for understanding the importance 
of dispersal time across a landscape matrix in affecting metacommunity dynamics. They further high- 
light the importance of landscape and dispersal patterns for predicting the onset of synchrony between 
weakly-coupled populations. 

Keywords Synchronization • Phase model ■ Dispersal • Time delay 
1 Introduction 

Predicting the onset and maintenance of synchronous fluctuations of abundance among populations has 
led to important progress in the understanding of population persistence, species coexistence, and response 
to environmental change. However, most theoretical studies of synchrony have assumed instantaneous and 
passive movement of individuals among discrete habitats. In natural systems, discrete populations are 
typically embedded within a landscape matrix that constrains patterns and duration of movement between 
populations. It is thus important to predict the importance of such a landscape matrix on the emergence 
and stability of synchronous fluctuations of abundance between distant populations. Here we incorporate 
delayed dispersal as a means to study landscape and dispersal patterns between discrete predator-prey 
communities. 

Synchronization is generally understood as a phenomenon where through interactions with one an- 
other, the timings of fluctuations in components of a dynamical system lock to a single pattern of variation 
in time. This concept arises in a variety of systems in physics, chemistry, social sciences, and biology: from 
crickets chirping in synchrony, to fireflies able to all flash in unison (Blasius and Tonjes 2007). Spatial 
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synchronization in the context of population ecology refers to coincident changes in population charac- 
teristics such as abundance, reproduction, mortality, or mean size or age distribution, in geographically 
separated populations (Liebhold et al. 2004). Ecological examples of synchronization include the cycling 
of hare-lynx populations, where fluctuations of groups across Canada were found to be locked on the same 
tight ten- year cycle (Blasius and Tonjes 2007). Synchronized population dynamics have been observed in 
natiiral systems across a variety of taxa and spatial scales, from fungal plant pathogens across 0.5-3m, 
insect herbivores across 1-lOOOkm, to birds across 5-2000km (Liebhold et al. 2004). For our purpose we 
study the synchrony of coupled oscillators, where individual populations undergo intrinsic oscillations 
in abundance with some natural frequency, and adjust their oscillations due to coupling (dispersal) be- 
tween oscillators. This formalism has been used to predict the onset of synchrony in simple ecological 
systems (Goldwyn and Hastings 2008). Networks of time-delayed and weakly connected oscillators have 
also been studied in physics, engineering, and neurophysiology (Schuster and Wagner 1989, Dhamala ct 
al. 2004, Campbell and Kobelevskiy 2012). However, no study has elucidated the role of delayed dispersal 
in coupled ecological systems. 

Causes of spatial synchrony have been studied extensively across systems and scales. Coupled oscillator 
theory predicts that oscillators can be synchronized through direct coupling, or by external forcing. 
Equivalently in ecology, there are two main mechanisms that produce synchrony: (a) density-dependent 
direct interactions, namely dispersal between populations, and (6) what is known as the Moran eflFect 
(Moran 1953), which is the entrainment of systems with similar density-dependent dynamical structures 
by correlated density- independent external factors, such as climate or weather (Liebhold et al. 2004). 
Because most natural systems arc connected through both dispersal and external stochastic fluctuations, 
and because both mechanisms produce similar patterns of synchrony it is often difflcult to determine 
which factor or combination of factors underlies observed patterns of spatial synchrony. 

Theory predicts that populations oscillating due to the same or slightly different density-dcipcndcint 
processes (linear or nonlinear) can be synchronized by dispersal of just a very small number of individ- 
uals per generation (Liebhold et al. 2004). However, if density-dependent processes are so different that 
populations oscillate with very different frequencies, then synchronization may not be possible through 
dispersal (Ranta et al. 1998). Besides dynamical stability of synchronization fluctuations, the rate of syn- 
chronization is another constraint on the ability of dispersal to explain synchrony in natural populations: 
dispersal-induced synchrony must converge rapidly to be observed (Goldwyn and Hastings 2008). Other- 
wise, populations will be kept in transient non-synchronous states by environmental perturbations. While 
strong coupling (i.e., a large number of species disperse per generation) can always increase convergence 
rates to synchrony in general, in a two predator-prey patch system connected only weakly by dispersal, 
the required property for fast convergence to synchrony is the separation on time scales between predator 
and prey dynamics (Goldwyn and Hastings 2008). This property is characteristic of relaxation oscillators 
where a large part of the cycle is spent at low population densities. In contrast, sinusoidal oscillators take 
a very long time to synchronize by weak dispersal. Our study examines some mechanisms underlying this 
result and uses its robustness to delayed dispersal to further predict characteristics of ecological systems 
that favour weak dispersal-induced synchronization. 

Metapopulation theory and the study of fragmentation have emphasized discrete boundaries and 
movement between populations (Levins 1969, Hanski 1999). The appeal of metapopulation theory lies in 
its simplicity associated with the assumption of passive and instantaneous population movement between 
patches. As a result, model predictions applied to natural systems usually ignore the complex landscape 
matrix that provides the context for discrete populations (Brady et al. 2011, Turner et al. 2001): moun- 
tain peaks embedded in valleys, forest patches within meadows and seagrass and between coral reefs. 
Dispersing propagules and migrating individuals can spend significant time and show non random move- 
ment within such a landscape matrix depending on the nature of the movement and on distance between 
populations. Current metapopulation theory assumes that the spatial structure of metapopulations can 
be captured by a per capita dispersal rate alone, rather than dispersal time. 

To better understand the mechanisms leading to synchronization in natural populations, we study the 
role of weak and delayed dispersal in driving the ecological conditions for the existence and stability of 
synchronous states, and for their resilience (convergence time). We illustrate how time delayed dispersal 
can be used to implement the particulars of individual movement (through active or passive dispersal) 
in the space betwccm patches. While discrete patch models are useful for studying synchronization of 
spatially distributed systems from a dynamical systems point of view, they are associated with the 
assumption of instantaneous dispersal where individuals disperse with no transmission delay. By including 
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a time delay that retains more detailed information about the dispersal process, we extend the patch 
dynamics model to integrate properties of the landscape matrix that contains discrete populations. We 
next formulate an ecological model as a two-patch Rosenzweig-MacArthur predator prey model coupled 
through prey dispersal with a time delay. We reduce our model to a phase model using weakly-coupled 
network theory (Appendix A). Our results show that synchronization through weak dispersal is highly 
sensitive to dispersal time. We find that with slight variation in the delay, the overall influence of weak 
dispersal on the rate of convergence to synchrony can vary greatly. This means that dispersal time between 
spatially discrete systems is a key variable for assessing the role of weak dispersal as a cause of synchronous 
fluctuations in natural systems. We further extend recent results from Goldwyn and Hastings (2008) on 
synchronization between patches to predict the importance of such a matrix on the synchronization of 
weakly coupled predator-prey metacommunities. 



2 Time-delay approach to modeling dispersal 

Let us first articulate the relationship of patch dynamics with time delayed dispersal to other spatial 
dynamical frameworks. Durrett and Levin (1994) showed the importance of choosing an approach to 
modeling spatially distributed systems- mean field approaches, patch models, reaction-diffusion equations, 
or interacting particle systems- that captures the right level of spatial detail for understanding population 
dynamics. The time-delay model we will consider (see (j3.2l) ) can be viewed as a combination of two existing 
treatments of spatial dynamics distinguished in this work. The system neglecting the time delay (i.e., for 
s = in (j3.2|) ) corresponds to the patch model approach: it groups individuals of a species into discrete 
patches without any spatially-explicit structure informing where patches are with respect to one another. 
In these models each patch is treated with a mean-field approach where individuals interact equally with 
each other and the patch itself is described by ordinary differential equations, and patches are connected 
to all others equally by dispersal or migration at a constant rate (Durrett and Levin 1994). 

Patch models and mean-field approaches both lose the level of spatial detail required to account for 
the movement and local interactions of individuals both in the patch, and in dispersal between patches. 
The most common approach to modeling spatial systems working on a finer level of spatial detail are 
reaction-diffusion equations, in which infinitesimal individuals diffuse in space and undergo purely local 
interactions (reactions). The dispersal term in a patch model (say, d{hj{t) — hi{t))) can actually be derived 
from reaction-diffusion equations directly, through the discretization of continuous space, and we take a 
moment to illustrate this. 

A common starting point for modeling random dispersal are diffusion equations, which are continuum 
approximations of discrete random walk processes. Individuals moving this way distribute over one- 
dimensional space and time according to the solution of the simple diffusion equation for an uncorrelated 
random walk: 

where _D is a diffusion constant related to parameters like step size, the probability to make a turn at 
each step as opposed to not moving, and the time taken for a step or a pause. If we shift from continuous 
one-dimensional space to discrete space so that Ui is the density at point in space x^, z = 1, . . . , n, (|2.ip 
becomes (for n=2): 

^ = 7jhW-".W], i,J = l,2-i^3 (2.2) 

where I is the distance between the patches. As shown by Allen (1983), this is the result of simply 
discretizing the second derivative in (j2.ip and applying Neumann boundary conditions, which confine the 
individuals to the union of the two patches, allowing no movement in or out of them. Equation (|2.2[) is 
the patch model representation of dispersal: it distributes u{x,t) over points in space representing the 
patches, so that information about the movement process between the patches is lost. It resembles the 
dispersal term in our main model p.2p , but without the incorporation of a time delay s for movement in 
the space between patches on the landscape matrix. 

By treating space as continuous, our time-delay model avoids discretizing intcr-patch space, keeping 
closer to the reaction-diffusion equation to describe dispersal in the space between patches. The addition 
of continuous space between patches (see Figure [T]) as opposed to discretizing space allows for handling 
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the concept of diffusive dispersal closer to the more ecologically realistic formulation based on finite- 
speed random walks in the environment separating the patches, so that dispersers take some meaningful 
amount of time to cross from patch to patch. The model retains key features of actual movement as it 
makes a random walk between the patches through the value of the delay r. Such features are lost in the 
discretized version with no delay. Furthermore, our time-delay dispersal is equivalent to the simplest case 
possible of the process of individuals reaction-diffusing in space: it assumes a constant mean dispersal 
duration, and no birth or mortality during dispersal. Under these assumptions we can equate incoming and 
outgoing population fiuxes from the landscape matrix to the patches themselves, thereby incorporating the 
dispersal process into patch dynamics equations. The assumption of no birth or mortality during dispersal 
is not necessary for the analysis in Sections |3][S] however, it is crucial for our numerical computations (in 
particular for the computation of H using XPPAUT in Section [S]). 

In a system of patches connected by dispersal, the mean transmission time r can be calculated in 
several ways using specific details of the individual dispersal process which might be easier to measure in 
natural systems than the time delay t itself. We now demonstrate two such derivations. We start with 
a continuous space distribution of individuals dispersing in a one-dimensional connecting path between 
patches, according to some underlying assumptions about movement. We then uncover the net squared 
displacement of individuals in the path, which give the mean time an individual spends dispersing between 
two patches separated by a distance I. 

If we first assume individuals move passively at a constant speed v with two-way travel on one- 
dimensional paths connecting patches, we can easily see that the mean crossing time r in this case is just 
v/l. Suppose we consider the more general movement pattern integrating active dispersal modeled by the 
telegraph equation: 

T d^u du _ v^T d^u 
2W^~dt ~ ~d^' 

where v is the finite speed of individual movement, and T is the characteristic time of movement before 
a direction reversal (see explanation in Turchin (1998) or derivation in Othmcr and Dunbar (1988)). The 
mean squared displacement over time, |x^(i)|, of a group of individuals instantaneously released into the 
system at patch 1 is calculated as follows in Othmer and Dunbar (1988): 

\x'{t)\=v'T[t-^(l-e~''/^)]. (2.4) 

While the mean displacement \x{t)\ is not quite equal to the square root of the mean squared dis- 
placement, one can scale •\/|a;2(t)| by some factor C to obtain a better estimate for \x{t)\ than ^y\x^t)\ 
(Byers 2001). Then the mean displacement is then 



\x{t)\^C^v^T[t--{l-e~^^/n]. (2.5) 

Now set \x{t)\ in (|2.5p to I and solve for the corresponding value of time t. This defines the delay s 

which appears in p.2p as a function of I, T, and v. Note that the right hand side of (12. 5p is a strictly 
monotonically increasing function of t, so there is a unique solution. 



3 Ecological model 

Following Goldwyn and Hastings (2008), we use the Rosenzweig-MacArthur predator-prey model for 
single-patch dynamics: 

dH f H\ acPH 

= rH\ 1 — 



dt \ K) b + H' , . 

dP aPH „ ^'^■^> 

= mP. 

dt b + H 

Here H{t) and P{t) represent the size of the prey and predator populations respectively. Prey growth 
is modeled by logistic growth with intrinsic rate r and carrying capacity K, and the intake of prey by 
predators is modelled with a Holling Type II functional response specified by the parameters a and b and 
where the conversion ratio of loss of prey to increase in predators is c > 1. The predators have a linear 
mortality rate, with parameter m. 
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Connecting path (x) 




Fig. 1 Two discrete (homogeneous) patches connected by a one-dimensional (two-directional) path for dispersal. 

Now suppose we have two identical and spatially discrete predator-prey patches, each with dynamics 
following (j3.1l) . We consider the case where only prey disperse between patches with a fixed dispersing 
rate D. Naturally, the dispersing prey take some finite transmission time s to cross the space between 
the patches. When we further assume the absence of growth or mortality during dispersal, we can equate 
the incoming flux of individuals to a patch at time t with the outgoing flux of individuals at the opposite 
patch at a time t — s, before the present time. This allows us to incorporate the dispersal process into 
our four-equation model (where Pi and Hi represent the predator and prey populations in the ith patch) 
with just the term D{Hj(t — s) — Hi{t)) in the rate equation for the prey in the ith patch. 

The result is the following two-patch Rosenzweig-MacArthur predator-prey model coupled through 
prey dispersal with a time delay: 

dt V K J b + H, ^ ' ^ 



dP, aP,H, 



dt b + H, 



(3.2) 

mPj, i,j = 1,2; « ^ j. 



This is identical to the model studied by Goldwyn and Hastings (2008), except for the time delay of s 
time units that we include in the dispersal term, and that we allow only one of the species to disperse 
between patches. 

We nondimensionalize p.2p following Goldwyn and Hastings (2008) to reduce the number of param- 
eters in a way that will be useful later for adjusting the timescale separation of the predator and prey 
growth rates. Thus we work with the following resulting equivalent form: 



dh; 

dt e 

dpi _ Pihj 

dt ~ 1 + h. 



- fh,{l~ah,)-^^) +d{h,{t-T)-h{t)), 
^ip^, i,j = 1,2; i j^j, 



where we reuse the variable t for scaled time at, and make the following substitutions: 
hi = Hi/b, Pi = [ac/rb]Pi, r = as, a = b/K, ^ = m/a, e = a/r, d = D/a. 

The single-patch dynamics of the Rosenzweig- MacArthur model are well-understood (Goldwyn and 
Hastings 2008). There is a region of s-a-fi parameter space in which the dynamics produce a stable limit 
cycle, namely: a < 1 and n < This is the case we are interested in for studying synchronization 

dynamics. 

Furthermore, we are concerned with the part of e-a-fj, parameter space where the predator-prey 
oscillations are relaxation-like, since this leads to faster convergence to synchrony when the patches are 
weakly coupled, a requirement for weak dispersal to be the cause of synchronization in nature. Goldwyn 
and Hastings (2008) explains how the separation of timescales between the predator and prey leads the 
single-patch system to relaxation- like oscillations: in the nondimensionalized form p.3p . this is achieved 
by decreasing any one of e, a, or /i. Decreasing e increases the rate of intrinsic prey growth relative to 
that of the predator; decreasing a increases the carrying capacity for prey, enhancing the size and time 
between prey outbreaks; decreasing fi means slower predator mortality, which leads to longer times for 
the predator population to decrease to the point where the prey population spikes, so that these spikes 
happen less frequently (Goldwyn and Hastings 2008). In effect, the prey populations spend more time at 
low numbers and then grow more rapidly at spikes when e, a, and jj, are small. 
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4 Phase model 

The purpose of this section is to reduce our ecological model to a phase model more practical to analyze. 
We use and refer to concepts of weakly connected network theory which are outlined in Appendix A. 

Suppose we choose parameters s, a, and /x for our model (j3.3l) so that the identical uncoupled sub- 
systems (i.e., the two patches with 6 = 0) oscillate on an exponentially orbitally stable limit cycle 7(i), 
which means that solutions starting close enough to the limit cycle approach it exponentially fast in the 
limit as time goes to infinity. The coupling in p.3p is symmetric, so the system can be written in the 
general form for weakly connected networks as follows: 

Xiit) = F(Xi(<)) + 6WiXi{t),X2it - r)) 
Mt) = F{X2{t)) + 6W{X2{t),Xiit - r)), 

where 

X.i = {hi,pif, 

nX,)^(i(,..,l-a/..)-J^i^),:^-.P.)' (") 

W{X,{t),X,{t - r)) - {hj{t - r) - h,{t),Q) 



T 



Characterized by its shape, position and natural frequency each periodic solution 7 is an isolated 
closed curve through two-dimensional phase space, a distorted circle that can be parameterized with a 
phase variable 6 that increases by 2ti in one period (see Appendix A). 

When we consider the full system, the fact that the 7 are exponentially orbitally stable limit cycle 
attractors means that for small 5, coupling between the Xi only significantly affects the phase variables 
9i. We assume that the time delay r is of an order of magnitude of 2ti/ fi or less. Then, by the theorem 
in Appendix A we can reduce (j4.ip to the corresponding two-dimensional phase model: 



^ = + 5H{e2 -01-^) 



(4.2) 



where 'P = Qt mod 27r and 

H{x) = ^j^ ^{t)-W{^{t),-i{t + x/Q))dt, (4.3) 

with 7(i) solving the system 



Finally, we define the phase difference variable (j){t) as (f) = 9i — 62. This reduces the system (j4.2p 
further, to 

^ = G{cj)) := S[H{-cf, ~ nr) ~ H{cb - fir)] = SHdeiav{<l>) (4.5) 
where we just denote with fir. 



5 Phase model analysis 

So far we have reduced the full system (14.11) to a one-dimensional dynamical system with variable (f) that 
is directly related to synchrony: (14. 5p relates the parameters e, a, and ^ using 12 and H, as well as a 
time delay r, to the relative phase positions of the two oscillators in their limit cycle. Now suppose the 
dynamics (p{t) from any starting phase difference (pito) converge with time to some fixed phase difference 
(p* with G{(f)*) = 0. We refer to this phenomenon as phase locking, and say that the system is in-phase 
synchronized when phase-locked at </>* = 0, and asynchronized when phase locked at ip* ^ 0. 



Synchronization in ecological systems by weak dispersal coupling with time delay 



7 



We solve for H numerically using the numerical software XPPAUT (Ermentrout 2002). We first find 
the numerical solution of a single patch model with given parameters s, a, and fi. We then calculate 
fi, and the function ^{t) that solves (|4.4p (the iPRC, discussed in Appendix A). XPPAUT can then 
can approximate the corresponding coupling function H{x) (the calculation uses (|4.3p ). We obtain this 
function via its first eleven Fourier coefficients a„ and 6„, n = 0,...,10, which approximate 

10 

H{x) = y~^[a„ cos(na;) + 6„ sin(na;)]. 
Now by a simple calculation (as done in Kobelevskiy (2008)), we define 

10 

= — 2 sin(n(/)) [a„ sin(ni7T) + 6„ cos(ni7T)] , 

Tl = 

so that calling c„ — an sin(n]7T) + 6„ cos(ni7T), we have an expression for G{(j)) = 

10 

G{(f)) = -26 Cn sin(n0). (5.1) 

n=0 

Our approximation was not improved by using more than eleven Fourier coefficients. We can now integrate 
G((/>) in time to find the solution of the system for an arbitrary starting phase difference to investigate 
the effect of r on phase dynamics. 

6 Results 

6.1 Phase-locked states 

We begin by investigating the effects of a time delay for the parameter set e = 0.1, a = 0.3, and /i = 0.3 
in p.3p . which for an uncoupled patch (i.e., S — 0) leads to stable oscillatory dynamics with frequency 
f2 = 0.74199. We use S — 0.001 for weak coupling. We retrieve our coefficients a„ and fe„ to construct 
G{(f)) ~ ^ from (|5.ip . which has roots corresponding to steady states. Notice that G((/)) is T-periodic in 
T. Thus wc examine a range of r values from to T = 8.468 — 27r//2. This range is justified under the 
assumption that r is of the same order of magnitude as T or less. If C {(/)*) < at a steady-state </>*, then 
it is stable, while if G'{(j)*) > then the steady state is unstable. 

The number and stability of phase-locked solutions to the two-patch predator-prey system is strongly 
affected by delayed dispersal (Figure [3]). We start with r = 0, equivalent to instantaneous travel between 
patches, and see that the system has two stable steady state solutions, and tt, and two unstable steady- 
state solutions at intermediate phase difference values (Figure [2] and [3]) . As r is increased, and tt 
are always steady-state solutions because G{(f>) is 27r periodic and odd. We see that other steady-state 
solutions may exist, appearing or disappearing in a total of ten saddle- node or pitchfork bifurcations as 
T is varied from to T = 8.468 (Figure E]). 

Our bifurcation analysis reveals the complex response of phase dynamics to a time delay in dispersal. 
Super-critical pitchfork bifurcations occur at the steady state (j>* = when r = 2.172, and at the steady 
state (j)* = TT when t = 6.407. In the bifurcation aX t — 2.172, as can be seen in Figure [3Ka), the steady 
state at 0* = changes stability, and a branch of stable steady states bifurcates from 0* = with the 
new stable branch existing for r < 2.172 where the 0* = solution is unstable. There is thus a stable 
steady state at 0* = on one side of the bifurcation, and a stable steady state close to on the other 
side of the bifurcation. Thus at this bifurcation the 0* = stable steady-state is carried continuously 
away from (or the reverse). The behaviour at the other supercritical pitchfork bifurcation at (j)* — t: 
when T — 6.407 is similar, except we see that there are actually two branches that bifurcate from <j)* — t: 
in symmetric fashion. This is characteristic of all supercritical pitchfork bifurcations (and also occurs in 
the bifurcation at 0* = 0, r = 2.172 on noting that = and (p = 2tt represent the same phase). 

Sub-critical pitchfork bifurcations occur at 0* = when r = 0.0228, 4.381 and 6.223 and at 0* = tt 
when T = 0.1462, 1.99 and 4.257. The bifurcations at </>* = when r = 4.381 and 6.223 are clearly visible 
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pi/2 pi 3pi/2 

Fig. 2 e = 0.1, a = 0.3, = 0.3, and <5 = 0.001. G(<j>) for a few different values of r. 




012345678 0,02 04 0.06 08 01 12 14 16 18 



(a) (b) 

Fig. 3 e = 0.1, a = 0.3, = 0.3, and 5 = 0.001. r-Bifurcation diagram: (a) showing a full r range, (b) showing values for 
T close to 0. Stable steady-states are black; unstable steady-states are blue (appearing as grey in monochrome). 

in Figure[3ja), while the bifurcation at (/)* = 0, r = 0.0228 can be seen in Figure[3jb). These are the most 
dynamically disruptive type of pitchfork bifurcation, since a stable steady state exists only on one side of 
these bifurcation points. A system in nature passing through such a bifurcation would 'jump' abruptly 
from one steady-state to another (outside the region of the bifurcation point). 

Saddle-node bifurcations correspond to two steady states (one stable and one unstable) that collide 
and annihilate so that a stable steady state either appears or disappears. There are two such bifurcations 
in our system; the one at t = 0.0124 is clearly visible in Figure [3jb) occurring close to subcritical 
pitchfork bifurcation at t = 0.0228. The second saddle-node bifurcation occurs at r = 4.247, close to the 
subcritical pitchfork bifurcation at = tt when r — 4.257. The separation of the t values of these last 
two bifurcations is very small and they are difficult to distinguish in Figure [3Ja), but the form of the 
bifurcation branches is similar to that of shown in Figure [31[b) for the other saddle- node bifurcation. 
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Table 1 e = 0.1, a = 0.3, fi = 0.3, and 5 = 0.001. Stable steady states for varying t values. 



r Range Stable steady-states 



(0, 0.0124) 


0, TT 










(0.0124, 0.0228) 


0, 


< 


r 


< 


TT, TT, -(p* 


(0.0228, 0.1462) 


< 


r 


< 


TT, 


TT, -(p* 


(0.1462, 1.99) 


< 


r 


< 


TT, 




(1.99, 2.172) 


< 


r 


< 


TT, 


TT, -(p* 


(2.172, 4.247) 


0, TT 










(4.247, 4.257) 


0, 


< 




< 


TT, TT, -0* 


(4.257, 4.381) 


0, 


< 


r 


< 


TT, 


(4.381, 6.223) 


< 


r 


< 




-r 


(6.223, 6.407) 


0, 


< 




< 


TT, -(j)* 


(6.407, 8.468) 


0, TT 











The fact that there are ten bifurcations across the range of r where the stable steady states change 
indicates that the quahtative dynamics of the weakly-connected system, and therefore its ability to 
synchronize, is remarkably sensitive to the value of the time delay r. Increasing r from to T can produce 
a full range (i.e., from to 27r) of stable phase differences (j>* (Table [1]). This enriches our understanding 
and ability to predict the dynamics of the system beyond the r = case with only two stable steady 
states (0 and tt). When r > it is still possible to have and/or tt as a stable steady state, but other 
asynchronous stable steady states are available as well over some values of r. However, there are also 
ranges of r values where the steady states and tt become unstable, separated by stable asynchronous 
steady states. Overall, our bifurcation diagram (FigureE]) can be used to show that the fraction of r values 
(of the full range length T) where in-phase synchronization is possible (i.e., is a stable steady-state) is 
(0.0228 + (4.381 - 2.172) + (8.468 - 6.223))/T = 0.529. This means that approximately half of aU possible 
time delays (on the same order of magnitude as T or less) prevent in-phase synchronization given the 
parameter values we used in numerical simulations. 

In fact, the system loses its ability to synchronize in-phase (i.e., the steady state becomes unstable) 
even with a very small time delay. If t were increased from to just 0.0228 (a time delay of 0.0228/a in 
unsealed time), a population with in-phase synchronized dynamics would bifurcate to a non-zero stable 
steady state, decreasing the risk of simultaneous extinction in both patches. The phase difference would 
remain held away from in-phase synchronization as r is increased, all the way until t = 2.172. This 
sensitivity of the dynamics to small changes of t also means that if two patches are synchronized in- 
phase, their global (simultaneous) risk of extinction can be reduced by a change in the delay r of less 
than (8.468 — 6.223)/2 = 1.123 to bring the dynamics to a sub-critical pitchfork bifurcation where the 
system shifts to another non-zero steady state, or to a super-critical pitchfork bifurcation where the 
steady-state phase difference smoothly increases from zero. 



6.2 Rate of convergence to synchrony 

Previous studies by Goldwyn and Hastings (2008) have shown that in general, relaxation oscillators 
converge to steady-state dynamics faster than other types of oscillators (due to their iPRCs having 
greater maximum magnitudes). We find that within this category of relaxation oscillators, a time delay 
plays a significant role in further determining how strong a synchronizing force weak dispersal can be. 
We find that convergence strength to synchrony predicted by the model varies widely over the range of r 
values we consider, and in a way that indicates that the system is highly sensitive to small changes of r. 

The convergence rate to synchrony of a stable steady-state 0* is calculated as |G'((/)*)|, the absolute 
value of the first derivative with respect to of the function G(0), evaluated at </>*; or in other words, 
the absolute value of the slope of the graph of G(0) at the zero (/>* (the slope will be negative if (p* is 
a stable steady-state). The greater the convergence rate to synchrony, the shorter the time for transient 
dynamics around the steady-state; and the more likely it is to be observed in natural systems exposed to 
recurring perturbations away from their phase-locked state. 

Our analysis shows that weak, non-delayed dispersal is not in itself a strong mechanism of in-phase 
synchrony in terms of convergence strength, compared to delayed dispersal. This is revealed by the slow 
convergence to in-phase synchrony at t = compared to those systems with delayed dispersal (r > 0) 
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(a) (b) 

Fig. 4 £ = 0.1, a = 0.3, fj, = 0.3, and 5 = 0.001. Rate of convergence to stable steady-state (blue) or the next stable 
steady state 0* > (red) if is unstable, (a) shows the values for the full t range, (b) shows a zoom of (a) for values closer 
to T = 0. 

that have an in-phase stable equUibrium {e.g. t = 3.2 or t = 6.9; Figure 2]), More specifically, the rate 
of convergence to a stable (jf = 0, or to the steady state with the smallest phase difference, is faster 
with than without a time delay (Figure H]). This is because at r = the system is close to the pitchfork 
bifurcation observed at t = 0,0228 where loses its stability, and it is precisely at that bifurcation points 
that G'{(j)) = and that the convergence rate is minimum. Convergence rate is fastest at the centres of 
the intervals between consecutive bifurcations points: these are the r values where weak dispersal is able 
to best produce synchronized dynamics in nature. Because there is a bifurcation very close to r = 0, a 
slight increase in r from zero not only creates asynchronous phase-locking possibilities, but also increases 
the rate of convergence to these phase-locked states, and thus the predicted ability of weak dispersal to 
explain synchrony in nature. Thus, we have shown that for a specific set of parameters corresponding to 
oscillatory dynamics of a single predator-prey patch, even a very small time delay can produce dynamics 
that are markedly different from the predictions of the non-delayed model. We conducted a sensitivity 
analysis showing the robustness of this result to a broad range of predator and prey parameter values (e, 
a, and fi; see Appendix B, 

7 Discussion 

The importance of a landscape matrix in controlling movement of individuals and propagules between 
discrete populations has led to its integration into recent management and conservation practices, but 
most theoretical frameworks of patch dynamics are still based on instantaneous dispersal. Our time-delay 
approach to modeling dispersal in a patch model constitutes a simple way to incorporate details of indi- 
vidual movement through the landscape matrix within a patch dynamics framework. We show that our 
patch model with delayed dispersal provides a useful approximation of dispersal through continuous land- 
scapes. Our results show that dispersal time is a spatial process deserving of attention when considering 
synchronization dynamics. We more precisely show how delayed dispersal can greatly limit the ability 
of weak dispersal to synchronize predator-prey dynamics to in-phase dynamics, which put a system at 
highest risk of global extinction. Instead, delayed dispersal is predicted to promote asynchronous phase- 
locked dynamics across discrete populations. Our results have important implications for understanding 
causes of spatial synchrony in natural systems and for bridging gaps between simple patch models and 
landscape ecology, 

7,1 Knowing when to be discrete, and when to be continuous 

Patch dynamics with delayed dispersal can provide a very useful tool to approximate complex movement 
patterns across natural landscape matrices. This approach effectively embeds a patch model onto a 
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more concrete physical space resembling a fragmented landscape, where local patches are separated by 
significant space that individuals disperse across. The modeling technique for dispersal here is applicable 
to any such system in nature where patches arc not directly next to each other (in which case r would 
just be zero), even if there's only a small transmission time between patches: there is no limit in general 
to how sensitive a system can be to small transmission delays. 

The notion of time delays in dispersal can be extended to systems involving egg banks, seed banks, 
dia-pausing eggs, or dormant stage in general. It is usually thought that a dormant stage is a means of 
"dispersal through time" to increase survival through harsh conditions (Callaghan and Karlson 2002). 
However, many examples of dormant stages also involve spatial aggregation and thus movement between 
dormant and non-dormant stages (Rothhaupt 2000). Such situations where dormant propagules disperse 
through space are just an extreme example of the complex ways that dispersal through space and time 
are coupled and can explain important ecological phenomena such as persistence and coexistence through 
storage effects (Chesson 2000). Our work shows indeed that spatial aggregation comes hand in hand with 
temporal segregation, which provides a relevant and insightful approach to modeling population dynamics 
in spatial systems. 



7.2 Time-delayed dispersal and synchrony 

Identifying causes of synchrony remains a major challenge in many biological systems. Ecologists have 
provided both theoretical and empirical evidence for weak dispersal as an important cause of spatial 
synchrony between discrete populations (see Liebhold 2004 for a review), but the question remains: when 
can weak dispersal be the mechanism behind in-phasc; synchronization between natural populations that 
are typically exposed to frequent phase perturbations and embedded within complex landscapes? Goldwyn 
and Hastings (2008) showed that a requirement for weak dispersal to cause synchronization in a network 
of predator-prey patches in nature is that the predator-prey dynamics behave as a relaxation oscillator. 
Relaxation oscillations result from a separation of temporal scales between predator and prey dynamics 
and lead to fast convergence time to in-phase steady state synchrony. Within this region of parameter 
space, our study provides a new perspective on the ability of weak dispersal to cause synchronization 
across natural landscapes where dispersal between habitats is delayed. 

Our time-delayed dispersal approach gives a way of using information we may have about a system 
connected through weak dispersal (the distance between patches, spciod of the dispersers, etc.) to infer 
the importance of weak dispersal as an important factor behind synchronized dynamics. Our results show 
that the delay itself can have a number of implications for phase dynamics: systems may be carried away 
from in-phase synchronization and held at a non-zero phase-locking state. Over a range of delay values 
(r) the system may increase from two to four the number of stable phase-locked steady states, further 
reducing the parameter space resulting in in-phase synchronization. Our results also relate the value of 
T to the rate of convergence to in-phase synchrony and show how a time delay can greatly impede the 
ability of weak dispersal to explain the maintenance of in-phase synchrony in nature, and instead enhance 
phase-locked dynamics. The rate of convergence is important because in nature where systems frequently 
experience perturbations due to stochastic fluctuations in the environment, their steady state dynamics 
are unlikely to be observed unless they recover quickly relative to the frequency of perturbations. In 
agreement with more general bifurcation theories (Strogatz 2000), our results show that the convergence 
rate is minimum near r bifurcation thresholds, and is maximum between these bifurcation points. This 
is especially relevant because for a broad range of s-a-fj, values, in-phase synchrony bifurcates to phase 
locking at delay values very close to r = 0. Instantaneous dispersal dynamics is thus highly sensitive to 
the introduction of arbitrarily small dispersal time. 

Our results suggest that the combined knowledge of all stable phase-locked states and of their con- 
vergence strength is important to predict when weak dispersal can be a very important cause of syn- 
chronization for systems characterized by relaxation oscillations. An example of such a system are those 
involving insect outbreaks, which commonly produce large amplitude oscillations (Peltonen et al. 2002). 
In communities characterized by small amplitude and sinusoidal oscillations, such as in the hare- lynx 
system, weak dispersal is an unlikely cause of synchronization, and extra information about the trans- 
mission delay may not be useful. Such systems are more likely to be synchronized through the Moran 
effect (Goldwyn and Hastings 2008), or by strong dispersal as dispersal strength decreases convergence 
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time. (The convergence rate \G'{(j)*)\ in our model inherits the factor 5 from G{4>) - see (|4.5I) or (I5.ip - 
and so the convergence rate is proportional to 6 and the convergence time is proportional to 1/^.) 

7.3 Patch dynamics within landscape matrices 

Fragmented landscapes are typically formed of habitats, connecting corridors, and the overall landscape 
matrix (Turner et al. 2001). Dispersal time captures a number of landscape properties: spatial arrangement 
of habitats, the presence and effectiveness of corridors, and the resistance of the landscape matrix to 
movement (Koh et al 2010). The idea of movement time (time-delay) between habitat fragments as 
assumed in our study can inform on the role of habitat corridors in connecting wildlife refuges. The use 
of these corridors in fragmented habitats to assist the dispersal of endangered species and decrease the 
risk of regional extinction is still controversial (Brady et al. 2011), and could be resolved through the 
understanding of dispersal delays across landscapes. 

While the presence of a landscape matrix and of corridors involve a finite dispersal time between 
habitats, it can also affect demographic processes and behavior of dispersing individuals, as well as 
regulate the rate of dispersal and modulate the coupling strength between habitats (Koh et al 2010). 
The phase-model we adopted allows studying movement time within a patch dynamical framework, but 
it also assumes weak dispersal and no change in individual density or behavior during dispersal. Future 
work should investigate the importance of a transmission delay with stronger dispersal coupling using 
delay-differential equations. The integration of demographic processes such as mortality during dispersal 
would also improve the relevance of patch dynamical models for understanding synchrony across complex 
landscapes. We hope our modeling approach can integrate these more realistic assumptions that are key 
to conservation, while still contributing to a more general theory of patch dynamics. 
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Appendix A 

For ease of analysis in studying synchronization with the assumption of weak dispersal, we reduce the 
full two-patch model (j3.3p to a phase model using concepts from weakly connected network theory (see 
(Hoppensteadt and Izhikevich 1977)). Here we expand on the phase model reduction used in Section UJ 
to understand how the general principles of weakly connected network theory can be used on a general 
system viewed from a coupled oscillator perspective. A weakly connected network in general is any system 
of the form 

^=F,{X,)+SW,iXi,...,X„), ^ = l,...,n, (A.l) 

where each Xi{t) e M™ and (5 is a small parameter. (For now we do not consider a time delay in coupling.) 
We are interested in the case where each decoupled subsystem {6 — 0) 

^^F,{X.,), t = l,...,n (A.2) 

has an exponentially orbitally stable limit cycle attractor 7^ C M™. 

The limit cycle % being a periodic orbit of the system, has an associated period Ti and a natural 
frequency f2i — iir/Ti. As a closed curve through m-dimensional space, each 7^ can be parameterized 
with a phase variable 6i that increases by 27r in one period. With an arbitrary starting point pi G 7i, we 
define a mapping Pi : [0, 27r) — >• that takes 9i G [0, 27r) to the unique corresponding point on 7^ that 
is P,{e,) = Y,{9^/Q,) = r,(0,T,/2^), where Y,(t) g solves dY,{t)/dt = F,{Y,{t)) with r,(0) = p {Y,{t) 
is on the limit cycle 7^) (Hoppenseadt and Izhikevich 1977). 

The fact that the 7^ are exponentially orbitally stable limit cycle attractors means that for small 5, 
coupling between the Xi only significantly affects the phase variables 9i. For systems where oscillators 



have identical frequencies Di = ■ ■ ■ — Qn — f^, Malkin's Theorem says that solutions of the system (jA.l|) 
can be continuously mapped to solutions of the following canonical phase model defined on the n-torus 
T" = §1 X • • • X (Hoppensteadt and Izhikevich 1977): 

'^ = n + 5H,{ei-e,,...,e^-e,) + o{5^), i^\,...,n, (a.s) 

at 

where Hi are phase coupling functions, 



TJo 



H,{ei~e,,...,er,~e,)^^ I i,{t)-wMt + {ei~0,)/Q),...,^r.it + ier,~e,)/f2))dt (A.4) 



and 7i(i) solves the system 



^ = -DFMt)rzit) 
z(t)-7z'(<) = l. 
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The function 7i(i) is referred to as the infinitesimal phase response curve (iPRC) (Winfree 1980; 
Kuramoto 1984), measures the degree to which an arbitrarily short and infinitesimally small perturbation 
advances (positive valued) or slows (negative valued) the phase (Goldwyn and Hastings 2008). The iPRC 
can be thought of as a measure of the sensitivity of the oscillator to perturbations at each time t in 
[0,Ti): at times when the oscillator is affected most by perturbations (from dispersal), the iPRC has 
a greater maximum magnitude. Goldwyn and Hastings (2008) show that the greater the separation 
in predator-prey timescales {i.e., relaxation oscillator), the greater the sensitivity of the oscillator to 
perturbations, and the higher the maximum magnitude of the prey component in the iPRC (Table [5]). 
This explains why relaxation oscillators converge faster to synchrony than more sinusoidal oscillators. It 
also underlies the requirement for separated timescales between predator and prey for weak coupling to 
result in synchronization in natural systems (Goldwyn and Hastings 2008). 

Crucially, phase model reduction is also possible for weakly connected systems involving an explicit 
dispersal delay. The phase model for such systems is derived in (Hoppensteadt and Izhikevich 1977) and 
(Ermentrout 1994). For weak coupling and 0{1) delays (i.e., delays of the same order of magnitude as 
the period of oscillation or less) , the delays do not explicitly appear in the phase model, but rather result 
in an additional phase shift. The main theorem of weakly connected delayed systems is stated as follows, 
adapted from (Hoppensteadt and Izhikevich 1977) and (Izhikevich 2008): 

Theorem (Phase model with delayed coupling): Consider a weakly connected oscillatory network 
that has an explicit transmission delay, described by the system 

^^ = FdXi) + 5WdXi{t-r^,i),...,X,,{t-r],n); X, e E", i = l,...,n, (A.5) 

where the rjij are finite nonnegative real numbers, all 0{1). Suppose that each uncoupled system has an 
exponentially orbitally stable T -periodic limit cycle solution ji C M™. Then, the system can be reduced to 
the phase model 

do 

-^ = n + SH,iei~e,~<Pu,-.-,9n^e,-^u) + 0{S^), i = l,...,n, (A.6) 
dt 

where <Pji — firjji mod 2t: , and the Hi functions are defined by ( [^.^[ j. 
Appendix B 

We examine specifically how varying the parameters e, a and /i might affect the general results we showed 
for a specific set of parameters leading to oscillatory dynamics of a single predator-prey patch. Table [2] 
shows the parameter sets, as well as the maximum magnitude of the iPRC at these parameter values. 
As expected, decreasing a and ^ (we keep e = 0.1) increases the height of the iPRC. Figures [5] and |6] 
shows the bifurcation diagrams analogous to Figure [3] for a number of other parameter sets. Similarly to 
the parameter set used for our main analysis, varying r across the interval [0,T] changes the dynamics 
and rates of convergence to steady states as the dynamics bifurcate a number of times. The bifurcation 
diagrams close to r = in Figure 6 show that there is a range of parameter space where = is stable 
at T = 0, but very close to a bifurcation, so that with even a slight increase in r, the system can move to 
asynchrony with a much faster convergence rate. Investigating this slightly larger region of e~a~fi phase 
space indicates that our qualitative results for the first set of parameters are generally comparable to 
other parameter sets, and that decreasing a and fi (separating the predator-prey timescales and making 
the oscillator relaxation-like) mainly has an effect that shows up on the dynamics at t = 0. 
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Table 2 Parameter sets. 



Parameter e 
set 




a 








iPRC (adjoint) 
height 
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0. 


1 





.2 





.15 


z.zo X iU 
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0. 


1 





.2 





.3 


1.3x10^ 
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0. 


1 





.2 
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0. 


1 





.2 
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1.1 X 10'^ 
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0. 


1 
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.2 


Q , , T nlO 
O.O X iU 
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l.UxlO 
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0. 
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0. 
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.5 
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0. 
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.3 





. 1 


■7 n T r\14 
; .U X iU 
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0. 
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.3 
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l.UxlO 
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.3 





.3 
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12 


0. 
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.3 
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0. 
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.3 
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20 


14 


0. 
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0, 


.35 


0, 


.1 


4.5xl0i° 


15 


0. 


1 


0, 


.35 


0, 


.2 


1.8x10'' 


16 


0. 


1 


0, 


.35 


0, 


.3 
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17 


0. 


1 


0, 


.35 


0, 


.4 


12 


18 


0. 


1 


0, 


.4 


0, 


.05 


2.5x101'' 
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0. 
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.4 
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4.5xl0'5 
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0. 
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0, 


.25 


130 
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.35 
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(a) 1 (b) 5 (c) 9 (d) 14 (e) 18 (f) 22 




(g) 2 (h) 6 (i) 10 (j) 15 (k) 19 (1) 23 




(m) 3 (n) 7 (o) 11 (p) 16 (q) 20 (r) 24 























(s) 4 


(t) 8 


(u) 12 (v) 17 


(w) 21 


(x) 13 



Fig. 5 5 = 0.001; refer to Table[2]for parameter set values corresponding to numbers 1 through 24. t Bifurcation diagrams 
(full range of r from to T) through a region of e-a-fi parameter space. Stable steady-states are black; unstable steady-states 
are blue. 
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(a) 1 



(b) 5 



(c) 9 



(d) 14 



(e) 18 



(f) 22 



(g) 2 (h) 6 



(i) 10 



0) 15 (k) 19 



(1) 23 



(m) 3 (n) 7 (o) 11 (p) 16 (q) 20 (r) 24 



(s) 4 (t) 8 (u) 12 (v) 17 (w) 21 (x) 13 

Fig. 6 5 = 0.001; refer to Table|2]for parameter set values corresponding to numbers 1 through 24. r Bifurcation diagrams 
(r from to 0.3) through a region oi e-a-fi parameter space. Stable steady-states are black; unstable steady-states are blue. 



